arXiv: 1509.01179v2 [physics.chem-ph] 18 Nov 2015 


Holomorphic Hartree—Fock theory: an inherently mnltireference approach 

Hugh G. A. Burton and Alex J. W. ThorrQ 
University Chemical Laboratory, Lensfield Road, Cambridge CB2 lEW, United Kingdom 

(Dated: November 20, 2015) 

We investigate the existence of holomorphic Hartree-Fock solutions using a revised SCF algo¬ 
rithm. We use this algorithm to study the Hartree-Fock solutions for H 2 and and report the 
emergence of holomorphic solutions at points of symmetry breaking. Finally, we find these holo¬ 
morphic solutions for H 4 and use them as a basis for Non-Orthogonal Configuration Interaction at 
a range of rectangular geometries and show them to produce energies in good agreement with Full 
Configuration Interaction. 


INTRODUCTION 

The Hartree-Fock (HF) approximation provides the 
bedrock of modern wavefunction based methods in Quan¬ 
tum Chemistry. Representing the wavefunction as a sin¬ 
gle Slater Determinant constructed from the product of 
one electron spin orbitals, a Self-Consistent Field (SCF) 
method is used to minimise the energy with respect to 
variation in the basis orbital^ thus solving the non-linear 
Roothaan-Hall equations[l|, ^ ■ It is a mean-field method 
whereby the electron-electron repulsion experienced by 
one electron is the average of the repulsion from all other 
electrons, and hence its solutions do not describe electron 
correlation effects. This is readily observed when study¬ 
ing the Restricted HF solution to the ground state of H 2 
which does not produce the correct behaviour at dissocia¬ 
tion. On the simplest level such behaviour at dissociation 
may be corrected by using an unrestricted Hartree-Fock 
(UHF) method, thus allowing the spatial parts of the a 
orbitals to differ to that of the /3 orbitals. Beyond this, 
the solutions of the Hartree-Fock method may be used as 
a reference for further correlation methods. 

Configuration Interaction (Cl) methods provide the 
most powerful solution to the problem of electron corre¬ 
lation. These take a linear combination of excited state 
determinants — obtained by replacing occupied orbitals 
with virtual orbitals — to introduce a term depend¬ 
ing on the inter-electron distance into the wavefunction. 
The Full Cl method uses all possible replacement deter¬ 
minants and yields the exact solution to the Hamilto¬ 
nian in the Hilbert Space of the system. Inevitably this 
method is very computationally expensive, scaling facto- 
rially with the number of electrons N, Q and so the use 
of a limited number of excited determinants is common 
in methods such as CISD 0{N^) which uses just single 
and double excitations. 

Recently, one of us has studied the existence of multi¬ 
ple solutions to the SCF equations using a method known 
as Metadynamics Q. These solutions have been investi¬ 
gated by a number of other authors and it is thought 
that they may provide a good representation of the ex¬ 
cited states of molecules. These low-energy SCF solu¬ 
tions have already proved applicable for a Restricted Ac¬ 


tive Space Self-Consistent Field method, @ and a study 
on LiF and Ogll^ found that they show similar proper¬ 
ties to diabatic molecular states of these molecules. In 
particular they do not obey non-crossing rules and they 
maintain a similar electronic structure as the molecular 
geometry is varied. 
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FIG. 1: The lowest energy RHF and UHF solutions for Hj 
using a STO-3G basis set are plotted and the disappearance 
of the UHF solutions may be observed at around I.2A. This 
point is known as the Coulson-Fischer point. These states are 
used as a basis for NOGI and a discontinuity in the ground 
state NOCI energy at the Coulson-Fischer point results from 
the disappearance of the UHF. 

The SCF states are not in general orthogonal, but 
may still be used in a configuration interaction cal¬ 
culation, and these states provide a suitable basis for 
Non-Orthogonal Configuration Interaction (NOCI) cal¬ 
culations^ yielding adiabatic states which show similar 


^ If very many SCF solutions (of the order of the size of the Hilbert 
space) are used in NOCI then there may be linear dependencies. 
However, we propose to only use the low-lying SCF states, which 
are generally significantly fewer in number than the total size of 
the Hilbert space, so we not envisage this to be a problem in 
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properties to the states produced by more expensive 
Cl methodsfl^. In particular, these states reproduce 
avoided crossings [l^, conical intersections and we be¬ 
lieve them to maintain size-consistency because they are 
formed from a basis of size-consistent SCF solutions. By 
using this basis, NOCI differs from a truncated Cl, whose 
size-inconsistency derives from the truncation at a fixed 
excitation level for a single reference. 

For this method to be useful it is imperative that these 
SCF states exist across all geometries of the system, yet 
the coalescence and disappearance of states is widely re¬ 
ported for a number of molecules [III ll2|- Even for the 
modest H 2 , it is found that the lowest energy UHF state 
coalesces (Figure [T]) with the lowest energy RHF state at 
the Coulson-Fischer point and disappears at shorter bond 
lengths, preventing the use of these states as a basis for 
NOCI calculations. Whilst the total number of solutions 
to the SCF equations is unknown (and one particular 
study by Fukutome 0 predicts it to scale dramatically 
with system size), to be useful the number of solutions 
should remain constant at all geometries. Using a re¬ 
vised Holomorphic Hartree-Fock theory[l^ we have lo¬ 
cated holomorphic-UHF solutions (hUHF) for STO-3G 
H 2 at bond lengths shorter than the Coulson-Fischer 
point. These correspond exactly to the UHF solutions 
where these exist, hence providing a constant number of 
solutions at all geometries. We believe this to also be the 
case for larger molecules. 

To find these hUHF solutions, one must solve a holo- 
morphized Schrodinger equation whereby the construc¬ 
tion of the density matrix is adjusted by removing all 
complex conjugates, producing a non-Hermitian Hamil¬ 
tonian. Previous methods using non-Hermitian Hamil¬ 
tonians have been used to study metastable electronic 
states and Feshbach resonances iB-T^ whilst complex 
valued orbital coefficients have been used in conjunc¬ 
tion with a standard Hermitian matrix to explore sym¬ 
metry broken RHF states. [ 2 ^ Our proposed holomor¬ 
phic method combines a non-Hermitian Hamiltonian 
with complex orbital coefficients which are solved to find 
hUHF solutions. These solutions are a basis for NOCI 
which provide a good description of the system ground 
state. 

In this paper, we study the existence and properties 


of these hUHF states further and propose a holomorphic 
SCF algorithm as a procedure for finding such solutions. 
This algorithm is used to find holomorphic solutions for 
H 2 in a minimal STO-3G and 6-31G* basis set and these 
are then used as a basis for a NOCI calculation before 
being compared to the FCI solutions in the same basis 
set. Finally, we apply the method to H 4 ''' and H 4 in 
both a square and rectangular geometry and compare 
our results for H 4 to those of other correlated methods. 


HOLOMORPHIC SCF 


Under a conventional SCF algorithm, the initial step 
is to guess a set of trial coefficients for N orbitals (de¬ 
noted expressed in a basis of M functions (de¬ 

noted /r, z/,...). These are used to generate a density 
matrix which is in turn used to form a Fock matrix. 
This Fock matrix is diagonalised to produce a new set 
of orbital coefficients and the occupied orbitals are se¬ 
lected. Using this new set of coefficients, the density ma¬ 
trix may be regenerated and the process repeated until 
self-consistency is reached. The convergence of this pro¬ 
cess may be accelerated using Pulay’s Direct Inversion 
of the Iterative Subspace (DHS) [211 [22| extrapolation, 
which uses a linear combination of previously calculated 
Fock matrices to minimise a DHS error vector (defined 
as FP — PF). 

The holomorphic UHF solutions can be found by lo¬ 
cating the stationary points on the holomorphic energy 
surface. In practice this is just an adjustment to the 
Hartree-Fock energy functional, achieved by removing all 
complex conjugates. Defining the holomorphic density 
matrix as 


N 




we replace the normal density matrices in the conven¬ 
tional Hartree-Fock energy functional (see e.g. Ref. 0 for 
notation) with density matrices of this form to give a 
holomorphic energy functional, 


E = + 2 ^ - {ficr\vT)CP>^'' -f 

fiv fipcrr 


Throughout the holomorphic SCF procedure, we have orbital f, 
normalised the orbitals using a method which also re¬ 
moves the complex conjugate, and require that for each 


M 

i = xcf. 
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Removal of these complex conjugates causes the func¬ 
tional to become a polynomial in the holomorphic density 
mat rix 14l . Taking inspiration from the single variable 
case[l^, where the Fundamental Theorem of Algrebra 
guarantees that there exists a constant number of solu¬ 
tions at all geometries, we have investigated whether the 
number of solutions also remains constant in this holo¬ 
morphic generalisation of Hartree-Fock theory. 

When using the holomorphic energy functional, the 
holomorphic density and Fock matrices are no longer Her- 
mitian but instead complex symmetric allowing the Fock 
matrix eigenvalues to become complex and thus selec¬ 
tion of new occupied orbitals may no longer be made 
using the Aufbau principle. Instead we use the Max¬ 
imum Overlap Method developed by Gill et aZ.0, se¬ 
lecting the new orbitals as the ones with the greatest 
projection into the space spanned by the previous set. 
This also makes it much easier to follow excited state 
solutions across geometries. Being complex symmetric 
matrices, the left and right eigenvectors are the same so 
we arbitrarily choose to use the right eigenvectors . 

In summary, our proposed holomorphic SCF algorithm 
proceeds as follows: 

1. Begin with a (possibly complex) guess for coeffi¬ 
cients Cf. 

2. Form the one-particle holomorphic density matrix 

Gf 

3. Generate the a and /3 Fock matrices 

- “ParifJ-crWr)]. 

+Y.ar[P'''" -^Pcr{^iCr\vT)]. 

4. Generate new orbitals from the right eigenvectors 
of the Fock matrix. 

5. Form the holomorphized overlap matrix of old and 
new orbitals. 

0 _ old^y new0l^ 

This overlap matrix is used to select the new or¬ 
bitals through the Maximum Overlap Method. 

6 . Repeat until self-consistency is reached. 

We try this algorithm with a large number of differ¬ 
ent initial guess sampled randomly to find as many sta¬ 
tionary points as possible. Once convergence has been 
reached, the true energy of the holomorphic solutions 
may be found by orthonormalizing the orbitals and then 
using the normal energy functional; these states may 
then be used alongside standard RHF and UHF solu¬ 
tions as a basis for configuration interaction methods. 
This step appears to be a principal difference to the 
non-Hermitian SCF methods proposed by McCurdy and 
Head-Gordon [l^ [l^ whereby the “true” exterior scaled 
wavefunction is evaluated using a transformed variable. 


Since the known RHF and normal UHF states are also 
stationary points on the holomorphic potential energy 
surface, they are also found by the holomorphic SCF al¬ 
gorithm, and included in the CL 


COMPUTATIONAL DETAILS 


Throughout the computational work, the necessary in¬ 
tegrals were generated using a modified version of Q- 
Chem 4.3 2^ and the SCF algorithm was implemented 
with additional processing using SciPyjl^. All figures 
were plotted with matr)lotlibp6|. 


HOLOMORPHIC SOLUTIONS TO H, 


The use of excited SCF solutions as a basis for config¬ 
uration interaction calculations is not possible even for 
H 2 . Breakdown of the RHF solution to yield a symmetry 
broken UHF solution (where the a and /3 electrons have 
become isolated on separate atoms) gives the Coulson- 
Fischer llj point and an inconsistent number of states 
across all geometries. Using our holomorphic SCF al¬ 
gorithm we have found hUHF solutions with complex 
coefficients which appear to correspond to these UHF so¬ 
lutions beyond the Coulson-Fischer point. 

Under a STO-3G basis set (corresponding to one 
atomic orbital centred on each H atom) we used QChem 
to generate the Ug and (T„ symmetry orbitals which were 
then used as the basis set for the holomorphic SCF al¬ 
gorithm. We identified two RHF and two UHF solutions 
(each UHF solution is doubly degenerate) and these are 
shown on the left-hand graph in Figure [2l the lowest of 
each matches those found previously in Ref. M We 
note the existence of a higher doubly degenerate RHF 
state corresponding to H+ ... H~ and H“ ... at dis¬ 
sociation however we do not search for these as accurate 
results may still be achieved with the states plotted in 
Figure [2] At the Coulson-Fischer point, a state can be 
seen rising sharply out of the UHF solution: this is a pair 
of holomorphic solutions which have complex coefficients. 
Whilst the cusp in the real energy of the hUHF state 
may appear distressing, we have previously found that 
the holomorphic energy is smooth and continuous 1^ . In 
the Non-Orthogonal Cl (NOCI) method, the six hUHF 
solutions are used to form a Hamiltonian matrix which 
is then diagonalized to produce a new set of new energy 
eigenvalues. Significantly, the NOCI solutions under the 
STO-3G basis show a smooth ground state energy which 
corresponds perfectly with the lowest energy FCI state. 
However, this is to be expected as under this minimal 
basis set the Hilbert space is completely spanned by the 
hUHF and RHF solutions and thus the NOCI method 
becomes equivalent to the FCI. 
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FIG. 2: The hUHF, NOCI and FCI states are plotted for Hj in a STO-3G (left) and 6-31G* (right) basis. The holomorphic 
solutions can be rising out of the Goulson-Fischer point where the normal UHF solutions disappear. For the STO-3G basis, 
the NOCI states correspond to the FCI states since their basis spans the full Hilbert space. For 6-31G*, there is still a very 
good correspondence and the observed NOCI states are smooth. 


H 


H 



FIG. 3: In this model, the size of the rectangle is dictated by 
the parameter R whilst 6 quantifies the degree of asymmetry. 


With this is mind, we repeat the calculations using the 
larger 6-3IG* basis, corresponding to two atomic orbitals 
centred on each atom. The set of six hUHF solutions are 
still used as a basis for NOCI but now the Hilbert space 
contains 16 determinants. These results are shown on 
the right-hand graph of Figure [2l Despite the hUHF 
solutions spanning only a subspace of the Hilbert Space, 
the NOCI states remain smooth as predicted[i^ and each 
one corresponds closely to one of the FCI states. The 
symmetries of these states may be identified (from the 
bottom) as ^S+,3 s+,1E+,1E+. 


SYMMETRY BREAKING IN H 


2 - 1 - 


H 2 is not the only molecule which exhibits disappear¬ 
ing UHF solutions, in fact many diatomics such as F 2 also 
have a Coulson-Fischer point and symmetry breaking of 
states has been observed in yet more complex molecules. 
Recently, Cohen et al. 12| have shown that symmetry 


broken states exist in the ground RHF state of H^^ with 
a rectangular arrangement. Using our holomorphic SCF 
method, we used a STO-3C basis to investigate this RHF 
state and its symmetry broken UHF states and these are 
plotted in Figure 0] for a square geometry and a rectan¬ 
gle with aspect ratio 5:7. Under this model it is possi¬ 
ble to simulate the interaction of two HJ molecules with 
very large ratios representing essentially non-interacting 
molecules and the square representing a strongly inter¬ 
acting system. 

In our results, the ground RHF state can clearly be 
seen covering all geometries and giving a poor description 
of the system energy in the dissociation limit. For the 
rectangular geometry, this RHF state breaks down into 
two symmetry broken UHF states at different geometries. 
The lower energy UHF state corresponds to a state with 
the a electron localised over a pair of hydrogens sharing 
a short edge of the rectangle and the /3 electron on the 
opposite side, thus maximising the bonding interaction 
between the pairs of H atoms. In contrast, the higher 
energy UHF state corresponds to the state with the a 
and P electrons localised over a pair of hydrogens shar¬ 
ing a long edge of the rectangle. At bond lengths shorter 
than Coulson-Fischer points, the hUHF states emerge as 
predicted to provide a continuation to the normal UHF 
states. These show complex coefficients as expected and 
as seen for the hUHF states of H 2 . For the square geom¬ 
etry, both possible arrangements of the a and /3 electrons 
become equivalent and thus the states become degener¬ 
ate. This causes the single four-fold degenerate hUHF 
state shown emerging from the RHF state. 

Using these 5 states (I x RHF and 4 x hUHF) as 
a basis for non-orthogonal Cl produces the states plot¬ 
ted in blue. Significantly, these states are smooth and 
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FIG. 4: The hUHF and NOCI states are plotted along with the ground FCI state for using a STO-3G basis for a square 
(left) and rectangular (right) geometry with aspect ratio 5:7. Two hUHF states can be seen for the rectangular geometry, both 
doubly degenerate and corresponding to the expected UHF state beyond the Coulson-Fischer point. In the square geometry 
these two states become degenerate to give a single symmetry broken state with four-fold degeneracy. 
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FIG. 5: The lowest energy RHF (red and blue) solutions for 
H 4 in a STO-3G basis are plotted at R=1.7A and show clearly 
that the RHF cusp observed is a crossing of the two states, 
giving two-fold degeneracy at 90°. A NOCI with just these 
two states recovers a smooth lower (green) and upper (not 
plotted) state. The corresponding CCSD states (magenta and 
cyan) are also plotted and also show a crossing at this point. 


continuous across all bond-lengths. The ground Full Cl 
state is plotted in green and shows good correspondence 
to the lowest energy NOCI state, which itself is a dra¬ 
matic improvement on the RHF state. It is also worth 
noting that further symmetry breaking of the UHF solu¬ 
tions were identified at larger bond lengths and this may 
provide a interesting base for further study. 


HOLOMORPHIC STATES OF H4 


We now proceed to examine the SCF solutions of four 
interacting H atoms in the arrangement depicted by Fig¬ 
ure O This model has been used repeatedly to assess the 
performance of Coupled-Cluster and other methods ( 2 ^ 
l29| since at large 9 the model resembles well separated 
H 2 units, where the multireference character of the wave- 
function increases, whilst at 90° the four atoms form a 
square and the solutions become degenerate. The au¬ 
thors of such studies have repeatedly found what they 
refer to as a cusp in the RHF energy at 90° which then 
produces subsequent cusps in the Coupled-Cluster ener¬ 
gies generated using this RHF reference. 


We have replicated the calculations carried out by 
Scuseria et al. [28| for this system at R=1.70A but using 
a STO-3G basis instead of Dunning’s DZP basis set, and 
we plot the two lowest energy RHF states along with their 
corresponding CCSD energies and the FCI energy in Fig¬ 
ure [5] In the Scuseria study, only the lowest energy RHF 
state at each geometry has been used as the reference for 
the CCSD calculation, ignoring any higher energy solu¬ 
tions. It appears obvious to us that one cannot ignore 
these higher energy states since Figure [5] demonstrates 
that this cusp is by no means an error in the RHF but 
is a crossing of two solutions, found to be degenerate at 
90°. Furthermore, it is clear that the subsequent cusp 
observed in the CCSD calculation is a crossing of two 
states in an analogous manner to the RHF states. Using 
these two RHF states as a basis for our Non-Orthogonal 
Cl method returns a smooth curve without any cusp. We 
note that though the shape models the FCI curve well, 
the energy of this NOCI solution is signihcantly too high, 
indicating that the large contribution due to dynamical 
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FIG. 6 : The 18 hUHF states found for square H 4 in a STO-3G basis are plotted in red and used as a basis for generating the 
NOCI states (blue) for which the ground state corresponds well with the ground FCI state (green). The degeneracy of the 
hUHF states at R=1.7A from bottom to top are: 2 , 4, 8 , 2 and 2 . 


correlation is absent. Using our holomorphic SCF algo¬ 
rithm we identified a further 16 hUHF states of similar or 
lower energy to these RHF states and these are plotted 
along the 9 = 90° cross-section of the potential energy 
surface in Figure IHl Again we observe the emergence of a 
holomorphic state with complex coefficients and the same 
degeneracy at the coalescence point of a UHF state with 
an RHF state, conserving the total number of states. We 
apply the NOCI method to the basis set formed from the 
18 RHF and hUHF solutions to generate NOCI solutions 
plotted in blue. The ground state NOCI state shows good 
correspondence to the FCI ground state solution, despite 
the NOCI space only containing 18 determinants com¬ 
pared to the 36 determinants of the Hilbert space. As 
a control we have compared to the NOCI formed from 
the RHF and 17 random complex UHF states (see sup¬ 
plementary Figure 1) where we have found the minimum 
deviation from the FCI to be Further to this, 

we have also plotted cross-sections through the poten¬ 
tial surface between 70 - 110° at R=0.65A, 0.85A and 
I. 70 A using the same hUHF states and these are shown 


in Figure [T] These cross sections contain 12, 8 and 0 
complex holomorphic states at 90° respectively. The two 
lowest NOCI solutions calculated from this hUHF basis 
are smooth and correspond well to the two lowest FCI 
states, again showing no cusp at 90°. It is also interesting 
to note that as one moves away from 90° in the R=0.85A 
plot, we see the 4 fold degenerate UHF turn holomor¬ 
phic whilst in the R=0.65A we see the 2-fold degenerate 
UHF state also yield a holomorphic state. Overall, each 
hUHF solution coalesces with an RHF solution at some 
point on this potential energy surface. In future stud¬ 
ies we hope to map these coalescence points across the 
potential energy surface. 


CONCLUSION 

In this study we have demonstrated that using a 
matrix-driven holomorphic SCF algorithm, it is possi¬ 
ble to locate holomorphic-UHF solutions. These hUHF 
solutions form an extension to the conventional UHF so- 
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FIG. 7: Cross sections through the H 4 potential energy surface in STO-3G with changing 9 are plotted at R=0.65A, O.SSA 
and I. 70 A and the hUHF states are used to generate NOGI states which show good correspondence to the FCI energies. 
Significantly, all of the NOCI states plotted are continuous and smooth through 90° and mirror the FCI energies almost exactly. 


lutions where such states coalesce with an RHF state and 
disappear. Despite being located by finding the station¬ 
ary states on a holomorphic, complex-valued potential 
energy surface, the hUHF states have real-valued holo¬ 
morphic energy expectation values, E, although we are 
unsure why this is the case. Furthermore the hUHF 
states provide an excellent basis for Non-Orthogonal Cl, 
yielding states which provide an excellent representation 
of the ground FCI state at a much lower computational 
cost. Since the UHF solutions are size extensive, we be¬ 
lieve the NOCI states will also be and we hope to verify 
this in the future. 

Using this algorithm, we have been able to demon¬ 
strate that the cusps observed in the energy calculations 
on H 4 actually represents the crossing of two RHF states 
and the ground FCI state may again be well represented 
using the hUHF states as a basis for NOCI. In future 
studies, we hope to explore the nature of these holomor¬ 
phic solutions further and look to use the hUHF solu¬ 


tions as a reference for other correlation methods such as 
Coupled- Cluster. 
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FIG. 1: For H 4 (with 9 = 90°) in a STO-3G basis, we generate 17 random sets of UHF orbitals (equally likely to be completely 
real, imaginary or complex) expressed in the RHF basis, and use these and the RHF at every point along the symmetric stretch 
to perform NOCI. We do not perform any convergence of these states but instead keep their orbital coefficients constant. 







































